%% Manylabs robust predictive analysis
% Replicates Figure 5, Figure A-6, Figure A-9 
% Run ManyLabs_main_ecf_over_nu.m first

% Set max iteration number and tolerance
max_iter = 5000;
tol = 0.0001;

% Random seed
seed = 101;

%% Case 1: N=66, J=5 (All effects included, nu estimated by EB)
N = 66;
J = 5;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet1 = 'N=66, J=5';
estimates1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'C3:H68'); 
se1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'N3:S68');   

% Randomly assign baseline and validation studies
[hat_theta_1, hat_vartheta_1, rep_base_se2_1, rep_val_se2_1] = random_assignment(estimates1, se1, seed);

% Compute and plot ecf as a function of nu
[empirical_cov_list_robust_1, rcv_cell_1] = ecf_over_nu_robust(hat_theta_1, rep_base_se2_1, hat_vartheta_1, rep_val_se2_1, meta_nu_list);
cov_plot_1 = plot_cov_over_nu_robust(empirical_cov_list_1, empirical_cov_list_robust_1, meta_nu_list);
saveas(cov_plot_1, '../Results/Robust_nu_plot_N_66_J_5.png');

% Plot histograms for interval length ratios
rcv_full = vertcat(rcv_cell_1{:});
figure;
CI_len_ratio_plot_1 = histogram(rcv_full(:,5) ./ rcv_full(:,6), 'Normalization', 'probability');
set(gca, 'FontSize', 18);  
box off;
saveas(gcf, '../Results/CI_len_rat_plot_N_66_J_5.png');  

% Plot histograms for variance ratios
figure;
var_ratio_plot_1 = histogram(rcv_full(:,1) ./ rcv_full(:,2), 'Normalization', 'probability');
set(gca, 'FontSize', 18); 
box off;
saveas(gcf, '../Results/var_rat_plot_N_66_J_5.png'); 

%% Case 2: N=132, J=2 (All effects included, nu estimated by EB)
N = 132;
J = 2;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet2 = 'N=132, J=2';
estimates2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'C3:E134'); 
se2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'K3:M134'); 

% Randomly assign baseline and validation studies
[hat_theta_2, hat_vartheta_2, rep_base_se2_2, rep_val_se2_2] = random_assignment(estimates2, se2, seed);

% Compute and plot ecf as a function of nu
[empirical_cov_list_robust_2, rcv_cell_2] = ecf_over_nu_robust(hat_theta_2, rep_base_se2_2, hat_vartheta_2, rep_val_se2_2, meta_nu_list);
cov_plot_2 = plot_cov_over_nu_robust(empirical_cov_list_2, empirical_cov_list_robust_2, meta_nu_list);
saveas(cov_plot_2, '../Results/Robust_nu_plot_N_132_J_2.png');

%% Case 3: N=48, J=5 (three irreplicable studies excluded, nu estimated by EB)
N = 48;
J = 5;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet3 = 'N=48, J=5';
estimates3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'C3:H50'); 
se3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'N3:S50');   

% Randomly assign baseline and validation studies
[hat_theta_3, hat_vartheta_3, rep_base_se2_3, rep_val_se2_3] = random_assignment(estimates3, se3, seed);

% Compute and plot ecf as a function of nu
empirical_cov_list_robust_3 = ecf_over_nu_robust(hat_theta_3, rep_base_se2_3, hat_vartheta_3, rep_val_se2_3, meta_nu_list);
cov_plot_3 = plot_cov_over_nu_robust(empirical_cov_list_3, empirical_cov_list_robust_3, meta_nu_list);
saveas(cov_plot_3, '../Results/Robust_nu_plot_N_48_J_5.png');

%% Case 4: N=96, J=2 (three irreplicable studies excluded, nu estimated by EB)
N = 96;
J = 2;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet4 = 'N=96, J=2';
estimates4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'C3:E98'); 
se4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'K3:M98');   

% Randomly assign baseline and validation studies
[hat_theta_4, hat_vartheta_4, rep_base_se2_4, rep_val_se2_4] = random_assignment(estimates4, se4, seed);

% Compute and plot ecf as a function of nu
empirical_cov_list_robust_4 = ecf_over_nu_robust(hat_theta_4, rep_base_se2_4, hat_vartheta_4, rep_val_se2_4, meta_nu_list);
cov_plot_4 = plot_cov_over_nu_robust(empirical_cov_list_4, empirical_cov_list_robust_4, meta_nu_list);
saveas(cov_plot_4, '../Results/Robust_nu_plot_N_96_J_2.png');
